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T— I Abstract 

Clustering with fast algorithms large samples of high dimensional data is an important chal- 
lenge in computational statistics. Borrowing ideas from |MacQueen ( 19671 who introduced a 
sequential version of the fc-means algorithm, a new class of recursive stochastic gradient algo- 
rithms designed for the fc-medians loss criterion is proposed. By their recursive nature, these 
algorithms are very fast and are well adapted to deal with large samples of data that are allowed 
to arrive sequentially. It is proved that the stochastic gradient algorithm converges almost surely 
0^ to the set of stationary points of the underlying loss criterion. A particular attention is paid to the 

,__( averaged versions, which are known to have better performances, and a data-driven procedure 

that allows automatic selection of the value of the descent step is proposed. The performance of 
the averaged sequential estimator is compared on a simulation study, both in terms of compu- 
tation speed and accuracy of the estimations, with more classical partitioning techniques such 
r* as fc-means, trimmed fc-means and PAM (partitioning around medoids). Finally, this new on- 

line clustering technique is illustrated on determining television audience profiles with a sample 
of more than 5000 individual television audiences measured every minute over a period of 24 



hours. 
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1 Introduction 



Clustering with fast algorithms large samples of high dimensional data is an important challenge in 
computational statistics and machine learning, with applications in various domains such as image 
analysis, biology or computer vision. There is a vast literature on clustering techniques and recent 
discussions and reviews may be found in Jain et al. ( 1999 1 or Gan et al. ( 2007[ ). Moreover, as argued 
in |Bottou| (|2010[), the development of fast algorithms is even more crucial when the computation 



time is limited and the sample is potentially very large, since fast procedures will be able to deal 
with a larger number of observations and will finally provide better estimates than slower ones. 

We focus here on partitioning techniques which are able to deal with large samples of data, 
assuming the number k of clusters is fixed in advance. The most popular clustering methods are 
probably the non sequential (Forgy ( 1965[ )) and the sequential (MacQueen ( 19671) versions of the 
fc-means algorithms. They are very fast and only require 0{kn) operations, where n is the sample 
size. They aim at finding local minima of a quadratic criterion and the cluster centers are given 
by the barycenters of the elements belonging to each cluster. A major drawback of the A;-means 
algorithms is that they are based on mean values and, consequently, are very sensitive to outliers. 
Such atypical values, which may not be uncommon in large samples, can deteriorate significantly 
the performances of these algorithms, even if they only represent a small fraction of the data as 



explained in Garcia-Escudero et al. (20101 or Croux et al. (20071. The A;-medians approach is a 



first attempt to get more robust clustering algorithms; it was suggested by MacQueen] ( 1967 1 and 



developed by [Kaufman and Rousseeuw (1990 1. It consists in considering criteria based on least 
norms instead of least squared norms, so that the cluster centers are the spatial medians, also called 



geometric or Li-medians (see Small ( 1990 1), of the elements belonging to each cluster. Note that 



it has been proved in Laloe (2010 1 that under general assumptions, the minimum of the objective 
function is unique. Many algorithms have been proposed in the literature to find this minimum. 
The most popular one is certainly the PAM (partitioning around medoids) algorithm which has 
been developed by Kaufman and Rousseeuw ( 1990| ) in order to search for local minima among 
the elements of the sample. Its computation time is 0{k-n?) and as a consequence, it is not very 
well adapted for large sample sizes. Many strategies have been suggested in the literature to reduce 
the computation time of this algorithm. For example subsampling (see e.g the algorithm CLARA 
in [Kaufman and Rousseeuw] (|1990|l and the algorithm CLARANS in |Ng and Han| (|2002[)), local 



distances computation (Zhang and Couloigner (2005 1) or the use of weighted distances during the 



iteration steps (Park and Jun (2008 1), allow one to reduce significantly the computation time without 
deteriorating the accuracy of the estimated partition. 



2 



Trimmed fc-means (see |Garcfa-Escudero et al.| (|2008[ |2010|) and references therein) is also a 



popular modification of the fc-means algorithm that is more robust (see Garcia-Escudero and Go- 



daliza ( 1999 1) in the sense that it has a strictly positive breakdown point, which is not the case for 



the A;-medians. Note however that the breakdown point is a pessimistic indicator of robustness since 
it is based on the worst possible scenario. For a small fraction of outliers whose distance is mod- 
erate to the cluster centers, fc-medians remain still competitive compared to trimmed A;-means as 
seen in the simulation study. Furthermore, from a computational point of view, performing trimmed 
/c-means needs to sort the data and this step requires O(n^) operations, in the worst cases, at each 
iteration so that its execution time can get large when one has to deal with large samples. 



Borrowing ideas from MacQueen ( 1967 1 and Hartigan ( 1975 1 who have first proposed sequen- 



tial clustering algorithms and Cardot et al. (2011 1 who have studied the properties of stochastic 
gradient algorithms that can give efficient recursive estimators of the geometric median in high di- 
mensional spaces, we propose in this paper a recursive strategy that is able to estimate the cluster 
centers by minimizing a /c-medians type criterion. One of the main advantages of our approach, 
compared to previous ones, is that it can be computed in only 0{kn) operations so that it can deal 
with very large datasets and is more robust than the A;-means. Note also that by its recursive nature, 
another important feature is that it allows automatic update and does not need to store all the data. 
A key tuning parameter in our algorithm is the descent step value. We found empirically that rea- 
sonable values are given by the empirical Li loss function. We thus also consider an automatic two 
steps procedure in which one first runs the sequential version of the /c-means in order to approximate 
the value of the Li loss function and then run our stochastic fc-medians with an appropriate descent 
step. 

The paper is organized as follows. We first fix notations and present our algorithm. In the third 
Section, we state the almost sure consistency of the stochastic gradient A;-medians to a stationary 



point of the underlying objective function. The proof heavily relies on |Monnez| ( |2006| ). In Section 4, 
we compare on simulations the performance of our technique with the sequential fc-means, the PAM 
algorithm and the trimmed /c-means when the data are contaminated by a small fraction of outhers. 



We note that applying averaging techniques (see Polyak and Juditsky ( 1992 1) to our estimator, with 
a small number of different initializations points, is a very competitive approach even for moderate 
sample sizes with computation times that are much smaller. In Section 5, we illustrate our new 
clustering algorithm on a large sample, of about 5000 individuals, in order to determine profiles of 
television audience. A major difference with PAM is that our algorithm searches for a solution in 
all the space whereas PAM, and its refinements CLARA and CLARANS, only look for a solution 
among the elements of the sample. Consequently, approaches such as PAM are not adapted to deal 
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with temporal data presented in Section 5 since the data mainly consist of and 1 indicating that 
the television is switched on or switched off during each minute of the day. Proofs are gathered in 
the Appendix. 



2 The stochastic gradient /c-medians algorithm 
2.1 Context and definitions 

Let (17, A, P) be a probability space. Suppose we have a sequence of independent copies Zi, . . . , Zn 
of a random vector Z taking values in W^. The aim is to partition i7 into a finite number k of clusters 
rJi, . . . , Qk. Each cluster Q,i is represented by its center, which is an element of M.'^ denoted by 0*. 
From a population point of view, the /c-means and fc-medians algorithms aim at finding local minima 
of the function g mapping M.'^^ to R and defined as follows, for x = {x^, . . . ,x^y with for all i, 
x' G M'^, 

g{x) e( min m\Z-x'\\)), (1) 

\r=l,...,k J 

where $ is a real, positive, continuous and non-decreasing function and the norm ||.|| in R'^ takes 
account of the dimension d of the data, for z G W^, \z\\' = Z]j=i -^j- The particular case 
= ti^, leads to the classical A:-means algorithm, whereas <\){u) = \u\ leads to the k-medians. 

Before presenting our new recursive algorithm, let us introduce now some notations and recall 



the recursive A;-means algorithm developed by MacQueen ( 1967 1. Let us denote by Ij- the indicator 
function, 

k 

Ir{z;x) = ^{\\z-x^\\<\\z-x3\\}^ 

i=i 

which is equal to one when x^ is the nearest point to z, among the set of points x*, i = I, . . . ,k. 



The /c-means recursive algorithm proposed by MacQueen ( 1967 1 starts with k arbitrary groups, each 



containing only one point, Xl , . . . , . Then, at each iteration, the cluster centers are updated as 
follows, 

-^^+1 = — ^nlriZn', Xn) {X''^ — Zn) , (2) 

where for n > 2, = (1 + n^)^^ and = 1 + Yl^^=i ^riZf, Xi) is just the number of elements 
allocated to cluster r until iteration n — 1. For n = 1, let = ^. This also means that X^_^_i is 
simply the barycenter of the elements allocated to cluster r until iteration n, 

^"+1 = l + ElJlAZt,Xe) Ui+fllr{Ze;X,)Z,] . 



1=1 
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The interesting point is that this recursive algorithm is very fast and can be seen as a Robbins-Monro 
procedure. 



2.2 Stochastic gradient A;-medians algorithms 

Assuming Z has an absolutely continuous distribution, we have 

P(||Z-x^|| = ll^-x^lD = 0, for stay j aadx' ^x^. 

Then, the fc-medians approach relies on looking for minima, that may be local, of the function g 
which can also be written as follows, for any x such that x^ ^ x^ when i ^ j. 



gix) = ^E[IriZ;x) \\Z - x' 



(3) 



In order to get an explicit Robbins-Monro algorithm representation, it remains to exhibit the gradient 
of g. Let us write g in integral form. Denoting by / the density of the random variable Z, we have. 



9{x) = yZl Ir{z;x)\\z - x'\\f{z) dz. 

7^1 J^-'Mx-} 



For j = 1, . . . , d, it can be checked easily that 



dx'^ 



X, — Zj 



and since 

Ir{z]x) 

the partial derivatives satisfy. 



xy 


-^3 


\\z — 





/(z)</(z), forz/x^ 



dg 
dx: 



{x) = / Ir{z;x)j^^—^f{z)dz. 



z — X' 



We define, for x G 



Vrg{x) = E 



IriZ;x) 



Ix"- - Zl 



(4) 



We can now present our stochastic gradient A;-medians algorithm. Given a set of k distinct 
initialization points in W^, X^, ■ • • , X^, the set of k cluster centers is updated at each iteration as 
follows. For r = 1, . . . ,k, and n > 1, 



-^n+l — — (^n^riZn, Xn) 



(5) 
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E 



— Zn 



■J n 



Fn = Zi, . . . , Zn-i). The steps aj^, also called gains, are supposed to be J>i-measurable. 

We denote by Vg{x) = {Vig{x), . . . , Vkg{x))' the gradient of g and define ¥„. = {V^, . . . V^)'. 
Let An be the diagonal matrix of size dk x dk, 
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each aj^ being repeated d times. Then, the /c-medians algorithm can be written in a matrix way, 

Xn+l =Xn- AnVg{Xn) - A„K, (6) 

which is a classical stochastic gradient descent. 

2.3 Tuning the stochastic gradient fc-medians and its averaged version 

The behavior of algorithm ^ depends on the sequence of steps aj^, r G {1, . . . , fc} and the vector of 
initialization Xi . These two sets of tuning parameters play distinct roles and we mainly focus on the 
choice of the step values, noting that, as for the fc-means, the estimation results must be compared 
for different sets of initialization points in order to get a better estimation of the cluster centers. 
Assume we have a sample of n realizations Zi, Zn of Z and a set of initialization points of the 
algorithm, the selected estimate of the cluster centers is the one minimizing the following empirical 
risk. 



1 " 

R{Xn) = - (7) 

i=\ r=l 

Let us denote by = 1 + YTi^ UiZf, Xi) the number of updating steps for cluster r, until 
iteration n — 1, for r G {1, . . . , k}. A classical form of the descent steps can be given by 



(1 + CaTlr 



if/,(Z„;X„) =0, 

Otherwise, 



(8) 



where c^, and 1/2 < a < 1 control the gain. 

Adopting an asymptotic point of view, one could believe that a should be set to a = 1 with 
suitable constants Cq and c^, which are unknown in practice, in order to attain the optimal para- 
metric rates of convergence of Robbins Monro algorithms (see e.g. Duflo ( \991) , Th. 2.2.12). Our 
experimental results on simulated data have shown that the convergence of algorithm (|5]) with de- 
scent steps defined in ([8]l is then very sensitive to the values of the parameters and Cq which 
have to be chosen very carefully. A simulation study performed in the particular case A; = 1 by 



Cardot et al. (2010l showed that the direct approach could lead to inaccurate results and is nearly 



always less effective than the averaged algorithm presented below, even for well chosen descent 



step values. From an asymptotic point of view, it has been proved in Cardot et al. (2011 1 that the 
averaged stochastic gradient estimator of the geometric median, corresponding to A; = 1, is asymp- 
totically efficient under classical assumptions. Intuitively, when the algorithm is not too far from 
the solution, averaging allows to decrease substantially the variability of the initial algorithm which 
oscillates around the true solution and thus improves greatly its performances. 



Consequently, we prefer to introduce an averaging step (see for instance Polyak and Juditsky 



(1992) or Pelletier ( 2000| )), which does not slow down the algorithm and provides an estimator 
which is much more effective. Our averaged estimator of the cluster centers, which remains re- 
cursive, is defined as follows, for r G {!,..., k}, n > 1, and for the value X^^^^ obtained by 
combining (|5]) and ([8]), 



^n+l 



UrXl + XI 



n.+l 



+ 1 



if/,(Z„;X„) =0, 
otherwise. 



(9) 



with starting points X\ = X 



1 ; 



1, . . . ,k. Then standard choices (see e.g. Bottou (2010 1 and 



references therein) for a and Cq are a = 3/4 and Cq = 1, so that one only needs to select values for 



3 Almost sure convergence of the algorithm 

3.1 A convergence theorem 

The following theorem is the main theoretical result of this paper. It states that the recursive algo- 
rithm defined in (|6]l converges almost surely to the set of stationary points of the objective function 
defined in ([3]l, under the following assumptions. 

(HI) a) The random vector Z is absolutely continuous with respect to Lebesgue measure, 
b) Z is bounded: 3K > 0: \\Z\\ < K a.s. 
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c) 3C: \/x e R'^ such that ||x|| <K + 1,E 

(H2) a) Vn > l.min^a;; > 0. 

b) max^ sup„ a"^ < min(^, g^) a.s. 

c) Yl'^=i maxj. aj^ = oo a.s. 

d) sup„ < oo a.s. 

(H3) Eti En=i «f < oo a.s. 



WZ-xl, 



< C. 



(H3') EtiEr=iiEk<)'/r(^„;X. 



< oo. 



Theorem 1. Assume that Xi is absolutely continuous and that \\X'[\\ < K,for r = 1, . . . , A;. Then 
under Assumptions (Hla,c), (H2a,b), (H3) or(H3'), g{Xn) and 

k oo 

r=l n=l 

converge almost surely. 

Moreover, if the hypotheses (Hlb) and (H2c,d) are also fulfilled then 'Vg{Xn) and the distance 
between Xn and the set of stationary points of g converge almost surely to zero. 

A direct consequence of Theorem 1 is that if the set of stationary points of g is finite, then the 
sequence necessarily converges almost surely towards one of these stationary points because 

Xn+i — Xn converges almost surely towards zero. By Cesaro means arguments, the averaged 
sequence Xn also converges almost surely towards the same stationary point. 



3.2 Comments on the hypotheses 

Note first that if the data do not arrive online and Xi is chosen randomly among the sample units 
then Xi is absolutely continuous and ||X[|| < i^, for r = 1, . . . , under (Hla,b). Moreover, the 
absolute continuity of Z is a technical assumption that is required to get decomposition Q of the 
Li error. Proving the convergence in the presence of atoms would require to decompose this error 
in order to take into account the points which could have a non-null probability to be at the same 
distance. Such a study is clearly beyond the scope of the paper. Note however that in the simple case 
/c = 1, it has been established in|Cardot et al.|(pOTT]l that the stochastic algorithm for the functional 



median is convergent provided that the distribution, which can be a mixture of a continuous and a 
discrete distribution, does not charge the median. 

Hypothesis (Hlc) is a stronger version of a more classical hypothesis needed to get consistent 



estimators of the spatial median (see Chaudhuri (1996 1). As noted in Cardot et al. (2011 1, it is 
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closely related to small ball properties of Z and is fulfilled when 

P(||Z-a;|| < e) < Ke^ 

for a constant k that does not depend on x and e small enough. This implies in particular that 
hypothesis (Hlc) can be satisfied only when the dimension d of the data satisfies d> 2. 

Hypotheses (H2) and (H3) or (H3') deal with the stepsizes. Considering the general form of 
gains aj^ given in ([8]l, they are fulfilled when the sizes rir of all the clusters grow to infinity at the 
same rate and a e]l/2, 1]. 



4 A simulation study 

We first perform a simulation study to compare our recursive /^-medians algorithm with the follow- 
ing well known clustering algorithms : recursive version of the A;-means (function kmeans in 
trimmed /c-means (function t kmeans in the package tclust, with a trimming coefficient a 
set to default, a = 0.05) and PAM (function pam in the ® package cluster). Our <® codes are 
available on request. 

Comparisons are first made according to the value of the empirical Li error (|7]l which must be 
as small as possible. We note that the results of our averaged recursive procedure defined by ([5]l, ([8]l 
and Q are very stable when the value of the tuning parameter c-y is not too far from the minimum 
value of the Li error, with a = 3/4 and = 1. This leads us to propose, in Section 4.2 an 
automatic clustering algorithm which consists in first approximating the Li error with a recursive 
/c-means and then performing our recursive fc-medians with the selected value of c-y, denoted by c 
in the following. We have no mathematical justification for such an automatic choice of the tuning 
parameter c but it always worked well on all our simulated experiments. This important point of 
our algorithm deserves further investigations that are beyond the scope of the paper. Note however 
that this intuitive approach will certainly be ineffective when the dispersion is very different from 
one group to another. It would then be possible to consider refinements of the previous algorithm 
which would consist in considering different values of tuning parameter c for the different clusters. 
We only present here a few simulation experiments which highlight both the strengths and the 
drawbacks of our recursive A;-medians algorithm. 
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4.1 Simulation protocol 
Simulation 1 : a simple experiment in 

We first consider a very simple case and draw independent realizations of variable Z, 

Z = (l-e)(7riZi+7r2Z2+7r3Z3) + e5„ (10) 

which is a mixture, with weights vri = 7r2 = vr3 = 1/3, of three bivariate random Gaussian vectors 
Zi, Z2 and Z3 with mean vectors /ii = (—3, —3), fi2 = (3, —3) and /X3 = (4.5, —4.5) and covari- 

/2l\ /3l\ [2-1 

ance matrices yar(Zi) = ,Var{Z2) = and Far (Z3) = 

Point z = (— 14, 14) is an outlier and parameter e controls the level of the contamination. A sample 
of n = 450 reaUzations of Z is drawn in Figure [T] 



▲ 




Figure 1 : Simulation 1 . A sample of n = 450 realizations of Z. An outlier is located at position 
(-14,14). 
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10 20 30 40 50 

Time index 

Figure 2: Simulation 2. A sample of n = 36 realizations of Z with d = 50. The mean values ^i, 
^2 and /i3 of the three natural clusters are drawn in bold hnes. 

Simulation 2 : larger dimension with different correlation levels 

We also performed a simulation experiment, with a mixture of three Gaussian random variables 
as in ([TOjl, but in higher dimension spaces with correlation levels that vary from one cluster to 
another. Now, Zi, Z2 and Z3 are independent multivariate normal distributions in M*^, with means 
/xij = 2sin(27rj7((i-l)),/X2j = 2 sin(27r/3 + 27ri/(d- 1)), and /isj = 2 sin(47r/3 + 27rj7(d- 1)), 
for j = 1, . . . , d. The covariance functions Cov{Zij, Za) = l.5p^^ for j, ^ G 1, . . . ,d and i G 
{1, 3} are controlled by a correlation parameter p, with pi = 0.1, p2 = 0.5 and pa = 0.9. Note that 
this covariance structure corresponds to autoregressive processes of order one with autocorrelation 
p. As before, 6z = (4, . . . , 4) S plays the role of an outlying point. A sample of n = 36 
independent realizations of Z, without outliers, is drawn in Figure[2]for a dimension d = 50. 
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4.2 Li error and sensitivity to parameter c 



As noted in Bryant and Williamson (19781, comparing directly the distance of the estimates from 
the cluster centers ^i, fi2 and //3 may not be appropriate to evaluate a clustering method. Our 
comparison is thus first made in terms of the value of the empirical Li error Q which should be as 
small as possible. For all methods, we considered that there were k = 3 clusters. 




Figure 3: Simulation 1 with e = 0.05 and n = 250. Mean empirical Li error (over 50 replications) 
for the PAM algorithm (dashed line), the fe-means (c = 0) and the stochastic A;-medians (bold line), 
for c g]0,10]. 



We first study the simple case of Simulation 1. The empirical mean Li error of the PAM 
algorithm, the A;-means and the averaged fc-medians, for 50 replications of samples with sizes n = 
250 and a contamination level e = 0.05 is presented in Figure[3] The number of initialization points 
equals 10 for both the fc-means and the /c-medians. When the descent parameter c equals 0, the 
initialization point is given by the estimated centers by the /c-means, so that the empirical Li error 
corresponds in that case to the A;-means error, which is sightly above 2.31. We first note that this 
Li error is always larger, even if the contamination level is small, than the PAM and the A; -medians 
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errors, for c e]0, 10]. Secondly, it appears that for c G [0.5, 4], the /c-medians Li error is nearly 
constant and is clearly smaller than the Li error of the PAM algorithm. This means that, even if the 
sample size is relatively small (n = 250), the recursive /c-medians can perform well for values of c 
which are of the same order of the Li error. 




"1 r 

1 



Parameter c 



Figure 4: Simulation 2 with n = 500, d = 50, and e = 0.05. The mean empirical Li error (over 
50 replications) is represented for the PAM algorithm (dashed line), the MacQueen version of the 
fc-means (c = 0) and the recursive /c-medians estimator (bold line), for c g]0, 7]. 



We now consider 50 replications of samples drawn from the distribution described in simulation 
2, with n = 500, d = 50 and e = 0.05. The number of initialization points for the A;-means and 
the A;-medians is now equal to 25 and the empirical mean Li error is presented in Figure |4] We 
first note that the performances of the PAM algorithm clearly decrease with the dimension. The 
/c-means performs better even if there are 5% of outliers and if it is not designed to minimize an Li 
error criterion. This can be explained by the fact that PAM, as well as CLARA and CLARANS, 
look for a solution among the elements of the sample. Thus these approaches can hardly explore 
all the dimensions of the data when d is large and n is not large enough. On the other hand, the 
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/c-medians and the A;-means look for a solution in all M'^ and are not restricted to the observed data 
and thus provide better results in terms of Li error. As before, we can also remark that the minimum 
error, which is around 1.36, is attained for c in the interval [0.5, 3]. 



<D - 



q 

CD - 



in - 




Parameter c 

Figure 5: Simulation 2 with n = 1000, d = 200, e = 0.05, and Z multiplied by a factor 10. The 
mean Li loss function (over 50 replications) is represented for the PAM algorithm (dashed line), the 
MacQueen version of the /c-means (c = 0) and our recursive fc-medians estimator (bold line), for 
c G]0,40]. 

We finally present results from Simulation 2 in which we consider samples with size n = 1000, 
of variable lOZ, with d = 200. The contamination level is e = 0.05 and 50 initialization points 
were considered for the /c-means and A;-medians algorithms. Since Z has been multiphed by a 
factor 10, the minimum of the Li error is now around 13.6. We remark, as before, that because of 
the dimension of the data, d = 200, PAM is outperformed by the fc-means and the /c-medians even 
in the presence of a small fraction of outliers (e = 0.05). The minimum of the Li error for the 
/c-medians estimator is again very stable for c G [5, 25] with smaller values than the Li error of the 
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fc-means clustering. 

As a first conclusion, it appears that for large dimensions the A;-medians can give results which 
are much better than PAM in terms of empirical Li error. We can also note that the averaged 
recursive A;-medians is not very sensitive to the choice of parameter c provided its value is not too 
far from the minimum value of the Li error. Thus we only consider, in the following subsection. 



the data-driven version of our averaged algorithm described in Section 2.3 in which the value of 
c is chosen automatically, its value being the empirical Li error of the recursive A;-means. This 
data-driven A;-medians algorithm can be summarized as follows 

1. Run the k-means algorithm and get the estimated centers. 

2. Set c as the value of the Li error of the k-means, evaluated with formula 

3. Run the averaged k-medians defined by ([5]), and ([9]), with c computed in Step 2 and Ca = ^■ 

4.3 Classification Error Rate 

We now make comparisons in terms of classification error measured by the Classification Error Rate 
(CER) introduced by [Chipman and Tibshirani ( 2005 1 and defined as follows. For a given partition 



P of the sample, let Ip(i^j') be an indicator for whether partition P places observations i and i' in 
the same group. Consider a partition Q with the true class labels, the CER for partition P is defined 

as 

^ ^ i>i' 

The CER equals if the partitions P and Q agree perfectly whereas a high value indicates disagree- 
ment. Since FAM, the fc-means and our algorithm are not designed to detect outliers automatically, 
we only evaluate the CER on the non-outlying pairs of elements i and i' of the sample. 

We present in Figure[6j results for 500 replications of Simulation 1, with a sample size n = 500 
and no outliers (e = 0). We note, in this small dimension context with no contamination, that the Li 
errors are comparable. Nevertheless, in terms of CER, the PAM, the A;-means and the data-driven 
fc-medians algorithms have approximately the same performances. For the trimmed fc-means, the 
results are not as effective, since this algorithm automatically classifies 5% of the elements of the 
sample as outliers. 

We then consider the same experiment as before, the only difference being that there is now 
a fraction of e = 0.05 of outliers. The results are presented in Figure [7] The /c-means algorithm 
is clearly affected by the presence of outliers and both its Li error and its CER are now much 
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Figure 6: Simulation 1 with e = and n = 500. On the left, the Li empirical error. On the right, 
CER for the A;-means, PAM, the data-driven recursive fc-medians algorithm (kmed) and the trimmed 
A;-means (tkm). 
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Figure 7: Simulation 1 with e = 0.05 and n = 500. On the left, the Li empirical error. On the 
right, CER for the fc-means, PAM, the data-driven recursive A;-medians algorithm (kmed) and the 
trimmed A;-means (tkm). 
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larger than for the other algorithms. PAM and the recursive A;-medians have similar performances, 
even if PAM is slightly better. The trimmed A;-means now detects the outliers and also has good 
performances. If the contamination level increases to e = 0.1, as presented in Figure [8j then PAM 
and the trimmed /c-means (with a trimming coefficient a = 0.05) are strongly affected in terms of 
CER and do not recover the true groups. The /c-medians algorithm is less affected by this larger 
level of contamination. Its median CER is nearly unchanged, meaning that for at least 50 % of the 
samples, it gives a correct partition. 
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Figure 8: Simulation 1 with e = 0.1 and n = 1000. On the left, the Li empirical error. On the 
right, CER for the fc-means, PAM, the data-driven recursive A;-medians algorithm (kmed) and the 
trimmed fc-means (tkm). 

We now consider Simulation 2, with a dimension d = 50 and a fraction e = 0.05 of outUers. The 
Li empirical errors and the CER, for sample sizes n = 500, are given in Figure|9] It clearly appears 
that PAM has the largest Li errors and the trimmed /c-means and the data-driven A;-medians the 
smallest ones. Intermediate Li errors are obtained for the fc-means. In terms of CER, the partitions 
obtained by the /c-means and PAM are not effective and do not recover well the true partition in the 
majority of the samples. The trimmed fc-means and our algorithm always perform well and have 
similar low values in terms of CER. 
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L1 error 
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Figure 9: Simulation 2 with e = 0.05, n = 500 and d = 50. On the left, the Li empirical error. On 
the right, CER for the /c-means, PAM, the data-driven recursive /c-medians algorithm (kmed) and 
the trimmed fc-means (tkm). 
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4.4 Computation time 

The '9 codes of all the considered estimation procedures call C routines and provide the same 
output. Mean computation times, for 100 runs, various sample sizes and numbers of clusters are 
reported in Table [T] They are based on one initialization point. From a computational point of 
view, the recursive A;-means based on the MacQueen algorithm as well as the averaged stochastic 
fc-medians algorithm are always faster than the other algorithms and the gain increases as the sample 
size gets larger. For example, when k = 5 and n = 2000, the stochastic A;-medians is approximately 
30 times faster than the trimmed A;-means and 350 times faster than the PAM algorithm. The data- 
driven recursive /c-medians has a computation time of approximately the sum of the computation 
time of the recursive /c-means and the stochastic A;-medians. This also means that when the allocated 
computation time is fixed and the dataset is very large, the data-driven recursive A;-medians can deal 
with sample sizes that are 15 times larger than the trimmed A;-means and 175 times larger than the 
PAM algorithm. 

Table 1: Comparison of the mean computation time in seconds, for 100 runs, of the different esti- 
mators for various sample sizes and number of clusters k. The computation time are given for one 
initialization point. 







n=250 






n=500 






n=2000 




Estimator 


k=2 


k=4 


k=5 


k=2 


k=4 


k=5 


k=2 


k=4 


k=5 


fc-medians 


0.33 


0.35 


0.36 


0.45 


0.47 


0.48 


1.14 


1.25 


1.68 


PAM 


1.38 


3.34 


4.21 


5.46 


15.12 


20.90 


111 


302.00 


596.00 


Trimmed fc-means 


2.20 


5.45 


6.76 


5.32 


11.19 


13.48 


22.97 


42.72 


51.00 


MacQueen 


0.21 


0.29 


0.31 


0.25 


0.43 


0.50 


0.60 


1.38 


1.76 



When the sample size or the dimension increases, the computation time is even more critical. 
For instance, when d = 1440 and n = 5422 as in the example of Section[5} our data-driven recursive 
/c-medians estimation procedure is at least 500 times faster than the trimmed /c-means. It takes 5.5 
seconds for the sequential fe-means to converge and then about 3.0 seconds for the averaged k- 
medians, whereas it takes more than 5700 seconds for the trimmed /c-means. 

5 Determining television audience profiles with /c -medians 

The Mediametrie company provides every day the official estimations of television audience in 
France. Television consumption can be measured both in terms of how long people watch each chan- 
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nel and when they watch television. Mediametrie has a panel of about 9000 individuals equipped at 
home with sensors that are able to record and send the audience of the different television channels. 
Among this panel, a sample of around 7000 people is drawn every day and the television consump- 
tion of the people belonging to this sample is sent to Mediametrie at night, between 3 and 5 am. 
Online clustering techniques are then interesting to determine automatically, the number of clusters 
being fixed in advance, the main profiles of viewers and then relate these profiles to socio-economic 
variables. In these samples, Mediametrie has noted the presence of some atypical behaviors so that 
robust techniques may be helpful. 
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Figure 10: A sample of 5 observations of individual audience profiles measured every minute over 
a period of 24 hours. 



We are interested in building profiles of the evolution along time of the total audience for people 
who watched television at least one minute on the 6th September 2010. About 1600 people, among 
the initial sample whose size is around 7000, did not watch television at all this day, so that we 
finally get a sample of n = 5422 individual audiences, aggregated along all television channels 
and measured every minute over a period of 24 hours. An observation Zi is a vector belonging to 
[0, l]*^, with d = 1440, each component giving the fraction of time spent watching television during 
the corresponding minute of the day. A sample of 5 individual temporal profiles is drawn in Figure 
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10 Clustering techniques based on medoids and representative elements of the sample (e.g. PAM, 



CLARA and CLARANS) are not really interesting in this context since they will return centers of 
the form of the profiles drawn in Figure [10] which are, in great majority, constituted of and 1 and 
are consequently difficult to interpret. Furthermore, the dimension being very large, d = 1440, 
these algorithms which do not consider all the dimensions of the data, as seen in the simulation 
study, will lead to a minimum value of the empirical Li error Q that will be substantially larger 
than for the fc-means and our recursive A;-medians. Indeed, at the optimum, the value of the Li 
empirical error is 0.2455 for the fc-medians, 0.247 1 for the fc-means and 0.2692 for PAM. 

The cluster centers, estimated with our averaged algorithm for A; = 5, with a parameter value 



selected automatically, c = 0.2471, and 100 different starting points, are drawn in Figure 11 They 
have been ordered in decreasing order according to the sizes of the partitions and labelled CI. 1 to 
C1.5. Cluster 1 (Cl.l) is thus the largest cluster and it contains about 35% of the elements of the 
sample. It corresponds to individuals that do not watch television much during the day, with a 
cumulative audience of about 42 minutes. At the opposite, cluster 5, which represents about 12% of 
the sample, is associated to high audience rates during nearly all the day with a cumulative audience 
of about 592 minutes. Clusters 2, 3 and 4 correspond to intermediate consumption levels and can be 
distinguished according to whether the audience occurs during the evening or at night. For example 
Cluster 4, which represents 16% of the sample, is related to people watching television late at night, 
with a cumulative audience of about 310 minutes. 
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Figure 1 1 : Cluster centers for temporal television audience profiles measured every minute over a 
period of 24 hours. 
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Appendix : Proof of Theorem [T] 



The proof of Theorem [T]relies on the following light version of the main theorem in Monnez ( 2006 1, 



section 2.1. The proof of Theorem [T] consists in checking that all the conditions of the following 
theorem are satisfied. 

Theorem 2 (Monnez (2006)). Assuming 

(Ala) g is a non negative function; 

{Alh) There exists a constant L > such that, for all n > 1, 

g{Xn^i) - g{Xn) < {Xn+i - Xn,Vg{Xn)) + L \\Xn+i - X„||^ a.s.; 

{Ale) The sequence is almost surely bounded and Vg is continuous almost everywhere on the 
compact set containing 

(^2) There exists four sequences of random variables (Bn), {Cn),{Dn) and (En) in M"*" adapted 
to the sequence {Tn) such that a.s.: 

{A2a) II VA:E[K|J-J||' < BrMXn) + C„ and EZii^n + Cn) < oo; 

{A2b) nWAnVnf \Tn] < D^giXn) + En and Zn=liDn + ^n) < OO; 

(^3) sup„ aj^ < min(|, ^) a.s., X^^j^ max^ aj^ = oo a.s. and 

maxr aZ, 

sup < cxD a.s. 

n miiir 

then the distance of Xn to the set of stationary points of g converges almost surely to zero. 
Proof of Theorem [7] 

Let us now check that all the conditions in Theorem 2 are fulfilled in our context. 
Stepl: proof of 

Let A = Xn and B = Xn+i. Since Xn is absolutely continuous with respect to Lebesgue 
measure, X]r=i ^r{Z; A) = 1 a.s. and one gets 

" k 

^Ir{Z;A) m.in\\Z -B^\ 



g{B) = E 



min \\Z — B 

_ r 



E 



,r=l 
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and it comes 



which yields 



g{B)<Y,nir{Z-A) \\Z-B-\\\, 



r=l 



g{B) - g{A) < J^E [7,(Z; A) {\\Z -B^\\-\\Z-A^ 



r=l 



The application x ^ Hz — x*"!! is a continuous function whose gradient 



T 

VII r II X — Z 

r\\Z — X — 



llx'' — z\\ 

is also continuous for x^ ^ z. Then almost surely for d>2, there exists = A^ + ^^{B^ 
< /x"" < 1, such that 

- 11-^ -^ni = (5''-^^v^iiz-cni). 

Consequently for all d>2, 

k 

g{B) - g{A) < J^E A){B- - A\ V. ||Z - CHI)] , 

r=l 

SO that 

k 

g{B)-g{A) < J2^[Ir{Z;A){B'-A',\7r\\Z-C'\\-\7r\\Z-A'\\)] 

r=l 
k 

+ ^ E [Ir{Z; A){B' - A^ \7r\\Z- A^)] ^ (1) + (2) 



On the one hand 



and on the other hand 



(2) = Y.{B^- ^^ VrgiA)) = {B-A, Vg{A)), 



r=l 



(1) < \\B'' -A''\\E [\\Vr \\Z -C^W-VrWZ - A 

r=l 

hence since 

llVr-ll^-Cni - Vr-ll^-^nill = 



C -Z M -z 



||C^-Z|| \A^-Z\ 



< 2 



P'--Z|| ' 



one gets, with (Hlc) 



(1) < 2^ WB"- -A'-W IIC'' - A^IE 



\Z- A"- 



< 2C^ US'" - A^f = 2C \\B - I 



r=l 
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Consequently, we have 



g{B) - g{A) < {B - A, Vg{A)) + 2C \\B - Af . 



Step 2: Proof of the assertion: Vn > 1, for all r = l,...k, \\ \\ < K + 2 sup„ a!^ 

Let us prove by induction on n that for all n G N*, for all r = 1, . . . , /c, ||X^|| < K + 2sup^ a^. 
This inequahty is trivial for the case n = 1: || < K. Let n G N* such that ||-^^|| < K + 
2sup„ a^, Vr G {1, . . . , k}. Let r G {1,. .. ,k}. First we assume that ||-^^|| < K + al^. Then it 
comes 

||^;+i|| < \m\+alIr{Zn;Xr,)<m\+a'^<K + 2a'^. 
Now in the case when K + < \\X^\\ < K + 2 sup„ aj^, one gets 

\\X'J>K + al>\\Z4+al, 

and then 

Since for Ir{Zn\Xn) = 0, = it remains to deal with the unique index r such that 

Ir{Zn; Xn) = 1. In that case, 

By (Hlb) and from the inequalities aj^/ — ^n|| < 1 and ||.^n|| < K < \\X^\\, we have, 

1 n \ II vru I „r ll^nll n vru 

' - \\xr-z4) "^-'1 + ""K-z^ii - 'l^-'l ' 

which leads to 1 1 X^^-^ 1 1 < ii' + 2 sup„ and concludes the proof by induction. 
Step 3: Proof of (yllc) 



From the integral form 



J 



^riz;x)j^ ^f{z)dz, 



it is easy to see that ^ is a continuous function of x. 
Step 4: Proof of (yl2a) 
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The definition of implies that j;^] = and hence E[F„| J'n] = 0. 
Step 5: Proof of (^26) 



E 



K 



fc 



r=l 



Hence assuming (H3), one gets 



E 



oo 



.n=l 



< OO. 



In the case when (H3') holds instead of (H3), one has 

oo k 



E 



.n=l 



< [{<flr{Zn;Xn)\ < oo. 



n=l r=l 



Consequently, 



which concludes the proof. 



oo 

EE [Pn^n II' l^n 



n=l 



< OO a.s, 



□ 
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